Simple lattice model for biological gels 
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We construct a three-dimensional lattice model for biological gels in which straight lines of bonds 
correspond to filamentous semi-flexible polymers and lattice sites, which are exactly four-fold coor- 
dinated, to crosslinks. With only stretching central forces between nearest neighbors, this lattice is 
sub-isostatic with an extensive number of zero modes; but all of its elastic constants are nonzero, 
and its elastic response is afBne. Removal of bonds with probability 1 — p leads to a lattice with 
average coordination number less than four and a distribution of polymer lengths. When bending 
forces are added, the dihited lattice exhibits a rigidity threshold at p = pi, < 1 and crossover from 
bending-dominated nonaffine to stretching-dominated affine response between pb and p = 1. 

PACS numbers: 87.16. Ka, 62.20.dc, 61.43.-j, 05.70.Jk 



Biological gels [11-3 ^^^^ elastic networks, formed by 
crosslinked semiflexiblc polymers, that play a critical role 
in determining and controlling the mechanical properties 
of eukaryotic cells. Here we introduce and analyze prop- 
erties of a three-dimensional (3d) lattice model for these 
gels in which straight sequences of bonds correspond to 
polymers, and lattice sites, with a maximum coordina- 
tion number of four, correspond to crosslinks. 

Much of our intuition about filamentous networks 
comes from studies of two-dimensional (2d) Mikado mod- 
els in which straight lines, representing semi-flexible 
polymers of length L with stretching modulus ji and 
bending modulus k, arc laid down randomly on a plane 
and crosslinked at their crossing points. These studies 
show that there is a crossover from non-affine, bending 
dominated to affine, stretching dominated response as 
the number of crosslinks per polymer is increased 0- 
01. The kagome lattice [Fig. [T]i] with coordination num- 
ber z = 4 is a periodic version of the infinite L limit 
of the Mikado model, albeit with a monodisperse dis- 
tribution of lengths between neighboring crosslinks (seg- 
ment lengths). This lattice with nearest- neighbor springs 
only and no bending energy exhibits a nonvanishing shear 
modulus [1] and affine response even though it is just on 
the verge of mechanical instability: With z = 2d { where 
d is the spatial dimension) under periodic conditions, it 
is exactly isostatic 0, , but it has a number of zero 
modes that scales as its perimeter Q. It provides a rig- 
orous demonstration of the existence of a lattice with an 
L — !■ cx) affine limit such as seen in simulations on the ran- 
dom lattice 0, Q ■ A network of crosslinked semi-flexible 
polymers in 3d still has a maximum of only four neigh- 
bors per crosslink, and it is subisostatic in the absence of 
bending forces with a number of zero modes that scales 
with its volume. It is, therefore, not obvious that the 
affine, stretching-dominated shear-rigid limit found in 2d 
can exist in such networks even though models assuming 
affine response are in good agreement with experimental 
measurements [lH. Indeed, 3d computer-generated fila- 
mentous networks [13, show bending but not stretch- 




FIG. 1. (Color online) Steps leading from a stack of kagome 
lattices (a) to our final model lattice (d). 



ing dominated linear response. 

By stacking and connecting kagome lattices, we con- 
struct a 3d model lattice with exact four-fold coordina- 
tion that supports shear and compressional stress even 
in the absence of bending forces. Using analytical theory 
and numerical simulations, we study the elastic proper- 
ties of this lattice as a function of the unitless measure 
k = K/(/xa^), where a is a length scale, of the relative 
strength of bending compared to stretching forces. We 
demonstrate analytically that our undiluted lattice does 
exhibit affine, bending- independent elastic moduli of or- 
der ^/a^ with that for pure shear in the xy plane equal to 
Go = We use numerical simulations to study the 

elastic properties of our lattice when polymers are short- 
ened by cutting bonds with probability 1 — p- Figure [2] 
provides a phase diagram summarizing our results. For 
all Pb < pi stretching (bending) dominates response at 
large (small) k. Near the rigidity percolation threshold 
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FIG. 2. (Color online) Schematic phase diagram for the 3d 
kagome lattice. lfp< pb, the lattice is floppy. Above pt, there 
are three regions: critical (I), dense (III), and transition (II). 
Each of these is divided into a bending dominated regime 
(clear) at low k and stretching dominated one (dotted) at 
high k. Ap= {p ~ Pb). 



Pb, G ^ {p — Pp)f with / w 0.2 with amplitude propor- 
tional to Go for K 1 and to n/a^ for k <^1. Near the 
dense limit p = 1, G/Gq is well described by a critical- 
like scaling function ofr~ p)^ ~ kL^, where 
L = a{l — p)^^ is the average polymer length, in which 
G approaches {hi/a^)L'^ for r ^ 1, but other scaling vari- 
ables cannot be ruled out. There is a transition region 
(II in Fig. [2]) that interpolates smoothly between between 
the rigidity percolation (I in Fig. [2]) and the dense (III in 
Fig. [2]) regions. These results are consistent with those 
found by Broedersz, Sheinman, and MacKintosh in 
a phantom lattice model, which like ours has a maxi- 
mum coordination of 4 and K-indcpendent affine moduli 
at p = 1. When k = 0, our simulations support the 
expected existence of a first-order jump in G at p = 1. 

We construct our 3d model lattice as illustrated in 
Fig. [TJ We start with a conventional kagome lattice with 
lattice constant a lying in the x-y plane with the x-axis 
along one of the straight lines of bonds. The kagome lat- 
tice has 3 sites in its unit cell as indicated by the shaded 
triangles in Fig. [TJ Next, we generate a stack of kagome 
lattices by placing replicas of the original lattice in planes 
normal to the 2;-axis with the spacing between consecu- 
tive planes equal to a. Then we generate 4 replicas of 
this stack. The first replica remains in its place. The 
second replica is rotated by an angle — 71/2 about the 
y-axis and then translated by a(l/8, 0, 1/4). The third 
replica is rotated by an angle </> = 7r/3 about the x-axis 
and then translated by a(0, •\/3/8,0). The fourth replica 
is rotated by an angle = — 7r/3 about the cc-axis and 
then translated by a(3/16, \/3/8,0). Finally, new lattice 
sites (crosslinks) are introduced at the crossing points of 
polymers. The resulting structure is shown in Fig.[l](d). 
It is a periodic lattice with strict 4-fold coordination and 
a tetragonal unit cell in the form of a rectangular paral- 
lelepiped containing 54 sites. As in Mikado models, the 
segment lengths in our model vary (from a/2 to a/16). 



The elastic energy density E' of a filamentous network 
is the sum of a stretching contribution and a bend- 
ing contribution E\^, each of which is a function of the 
displacements Uq, of lattice crosslinks a from their rest 
positions x^. To harmonic order in elastic displacement, 
we can write these contributions for our model lattice as 
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where Au^^ = - u,3 and Axa/j x^ - X/j. e^a,i3) 
is the unit vector directed from site a to /3 in the ref- 
erence lattice. The summations run over all bonds and 
bond-pairs, respectively. K-ya/s = 2k/(|Axq,^| -I- |Ax^a|) 
is a segment-length dependent bending constant derived 
from the worm-like-chain model. In the following, we 
treat ^ and k as independent mechanical parameters; in 
real systems p, is dominated by entropic stretching and 
is a function of temperature and k [15| . 

Under imposed external strain, individual lattice sites 
undergo displacements Ua.i = VijXa.i + Sua^i for i = 
X, y, z, where rjij is the imposed macroscpic deformation. 
When the equilibrium value of Su^ is zero, each dis- 
placement follows the macroscopic strain, and response 
is ajjine; when it is nonzero, response is nonajjine. In 
both the affine and nonaffine cases, the elastic energy 
£ = E/{Na^) per crosslink per volume of our model 
obtains the form appropriate to tetragonal symmetry 
with 6 independent elastic constants: 
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wher 



, = ^{pij+riji) is the usual, linearized symmetric 
Lagrange strain tensor [l^ and N is the number of sites 
in the lattice. 

We consider first the undiluted model with k = 0. In 
this case, a straightforward symbolic solution for all 
in terms of rjij yields Sua = for all a, i.e., response 
is purely affine. This implies that all filaments remain 
straight under elastic distortion, and, as a result, the 
clastic energy is independent of k. Our model thus pro- 
vides a proof of principal of the existence of 3d central- 
force, subisostatic lattices with purely affine response. 
The elastic constants of the undiluted lattice in units of 
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The polymers forming biological gels have finite length 
and they are polydisperse. Furthermore, the topology of 
their networks is that of a random solid rather than of a 
crystalline solid with a well defined point group symme- 
try. To study the infiuence of randomness and network 
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FIG. 3. (Color online) (a) Shear modulus Q; (b) fraction of 
rigid conformations n; (c) non-affinity parameter F. 



connectivity on the elasticity of our model lattice, we 
dilute it by randomly removing bonds with a given prob- 
ability 1 — p. Then, we calculate its mechanical response 
numerically for a range of values of p, k and number of 
crosslink sites N = 545, where S = S^SySz is the num- 
ber of unit cells stacked Saj-times in the a;-direction, etc. 
We focus on the response to shear in the x-y-plane and 
set r]ij = ^{SixSjy + SiySjx) with imposed strain 7 — 0.01 
for all deformations. We generate 100 random conforma- 
tions, and for each we calculate all displacements 6ua by 
minimizing the total energy using a conjugate gradient 
method. For each k, p, and TV, we calculate G, the stan- 
dard measure of nonaffinity S IH F = ^ {Su^^Y , 
averaged over all configurations, and the fraction n of 
configurations with nonvanishing G. 

Figure [3] shows plots of these quantities for S = 72 
unit cells as a function of p for different k. The ratio 
G ~ G{p, k)/Go approaches 1 as p — > 1, as expected, and 
it does so more rapidly with increasing k such that re- 



sponse is nearly affine for p > 0.8 for all but the smallest 
K, indicating that response for physical dense networks 
can be nearly affine. For all k > 0, G vanishes at a rigidity 
percolation threshold PbiS). For k = 0, however, G van- 
ishes at a much larger threshold Pc{S). Both pb{S) and 
Pc{S) decrease as the system-size decreases, with 5' -> 00 
values of = 1 and pb ~ 0.602, the latter in good agree- 
ment with the Maxwell counting arguments of Ref. [18| . 
The fraction n reaches unity, its upper limit, for k > at 
p ~ 0.6, a value that changes little with S. For smaller p, 
it drops down to zero with steepness that increases with 
S and that approaches a unit-step function for S* — ^ 00. 
The values of n between and 1 below p ~ 0.6 are a fi- 
nite size effect. The corresponding conformations are the 
ones that lead to nonzero values of G in Fig. [3^ below 
p ^ 0.6. For K = 0, n approaches a unit step function at 
p = 1 for 5 — > 00. The nonaffinity parameter F, shown 
in Fig. [3j: has a peak at p ^ 0.6 for all n and vanishes 
as required for p = 1; for p just below 1, however, it 
increases substantially for small k. The location of the 
peak changes little with S. 

Based on these observations, we expect the following 
scenario for the infinite-size limit: For k = 0, G displays 
a discontinuous jump from zero to its affine value a,t p = 
Pc = 1 reminiscent of a first-order phase transition. For 
all K > 0, G undergoes a rigidity percolation transition 
at Pb ^ 0.6. The behavior for both k ~ Q and k > is in 
qualitative agreement with the 2d kagome lattice, where 
the shear modulus displays a first-order jump at p = 1 for 
k = {) and a rigidity percolation transition &X, pb ^ 0.6 for 
K > Guided by effective medium theory (EMT) 

for the 2d kagome lattice [31, we fit G near pb to the 
scaling form 

g{k)\^p\f, 



Q{p, a) 
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where Ap ~ p~pb and g{k) — cik/(c2 + k). As shown in 
Fig. H] data-collapse is obtained for pb = 0.602, ci = 0.08, 
C2 = 0.1 and / = 0.2, but we cannot rule out somewhat 
different functions g{k) and values of / as small as zero, 
which leaves open the possibility of a weak first-order 
transition. Interestingly, the exponent / = 0.2 is not far 
from the exponent /3 = 0.175 for the size of the perco- 
lating rigid cluster in simulations of the transition to a 
rigid but unstressed state in models of network glasses 
[20| . Near p ~ 1 , the kagome EMT suggests the scaling 
form 

g(p,K) = T-i(-i + \/rT7)', (5) 

as long as the scaling variable t = Ak/{1 — pY is greater 
than lO'^K. With A = 0.7, this form provides a excel- 
lent fit to our data for t > 10 Kj BjS shown in Fig. [5] 
The EMT scaling form crosses over from G ^ Gq for 
T > 1 to G = iGoT ~ (K/a4)(l - p)-2 {k/a^)L^ 
for r <C 1, and we expect our model network to show 
the same crossover for 5 — > 00, although further simula- 
tions with larger system-sizes and smaller values of k are 
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FIG. 4. (Color online) Scaling of Q near pb- 



kagome lattice, compression moduli vanish for all -0 > 
but contrary to the kagome lattice, shear moduli remain 
nonzero only up to a small critical value ^pc oi ip. Thus 
for ijj > ipc, shear moduli at p = 1 vanish with k as they 
do in the diamond 2l[ and computer generated network 
lattices [11 [H. 
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FIG. 5. (Color online) Scaling of Q near p = 1. The symbol- 
legend is the same as in Fig|4l 



needed to corroborate this conjecture which is consistent 
with the findings for the 3d phantom lattice 
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all, the EMT predictions for the 2d kagome lattice work 
remarkably well for our 3d kagome based lattice. The 
likely explanation is that bending forces are quite effec- 
tive at restoring rigidity, and in the end, the important 
thing is that both the 2d and 3d lattices have nonvanish- 
ing bulk and shear moduli at p = 1 and k = and both 
exhibit a first-order rigidity transition at p = 1 when 
K = 0. 

In our network, all filaments are straight. In real net- 
works of semi-flexible polymers, filaments are in general 
not straight, and as a result, like the 4-coordinatcd dia- 
mond lattice 2l| , they may not have nonvanishing shear 
and bulk moduli [l3| when k — Q even in the limit 
p = 1 (L oo). If unit cells in the kagome lattice 
are twisted through an angle ip so that filaments are no 
longer straight, the lattice no longer resists compression 
when K = 0, but it does resist shear These examples 
make it clear that network geometry can play a roles as 
important as coordination number in determining elastic 
response. Further research on lattices with different bend 
geometries is clearly of interest. We have begun studying 
a twisted version of our 3d lattice, and our preliminary 
results indicate that at p = 1 and k = as in the twisted 
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